---
title: "scenario analysis"
output: html_document
date: "2024-02-06"
---

## Read me ##

This code has 4 components

1) It calculates and plots the GHG emissions per meal for the different fuel sources
2) It calculates the total emission from the rural and urban areas annually 
3) It calculates and plots the GHG emission for each sub-division in Senegal
4) Takes raster output from Khavari et al 2023 and plots the time saved after switching to clean cooking fuel for each sub-division in Senegal


# **All the files needed to run this code is included in the folder ** #

# ** Please change the file path to run the analysis **#

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("/Users/yusufjameel/Dropbox/ProjectDrawdown/Senegal/SEN/GIS_data")
library(reticulate)
library(sf)
library(raster)
library(ggplot2)
library(viridis)
library(RColorBrewer)
library(dplyr)
library(rnaturalearth)
library(rnaturalearthdata)
library(rnaturalearthhires)
```

## Biomass

```{r biomass, echo=FALSE}

##parameters
Ubiomass = 16 #MJ/Kg
CO2_biomass = 112 #Kg/GJ
CH4_biomass = 0.8640 #Kg/GJ
N2O_biomass = 0.0039 #Kg/GJ
BC_biomass = 0.1075 #Kg/GJ
OC_biomass = 0.3080 #Kg/GJ

GWP_CH4_biomass = 25
GWP_N2O_biomass = 298
GWP_BC_biomass = 900
GWP_OC_biomass = -46

fRNB = 0.285

##
gamma_biomass = CO2_biomass*fRNB + CH4_biomass*GWP_CH4_biomass + N2O_biomass*GWP_N2O_biomass + BC_biomass*GWP_BC_biomass + OC_biomass*GWP_OC_biomass ## emissions unit is kg/GJ

##
efficiency_traditional = .12
efficiency_nautraldraft = .30
efficiency_forceddraft = .37

##
energyuse_traditional = (3.64/efficiency_traditional)*10^-3
energyuse_nautraldraft = (3.64/efficiency_nautraldraft)*10^-3
energyuse_forceddraft = (3.64/efficiency_forceddraft)*10^-3

##
emissions_traditional = energyuse_traditional*gamma_biomass
emissions_forceddraft = energyuse_forceddraft*gamma_biomass
emissions_nautraldraft = energyuse_nautraldraft*gamma_biomass

#
wt_emissions_traditional  = 1
wt_emissions_forceddraft = 0
wt_emissions_nautraldraft = 0

##
mean_emissions_biomass = wt_emissions_traditional*emissions_traditional + wt_emissions_forceddraft*emissions_forceddraft + wt_emissions_nautraldraft*emissions_nautraldraft

```

## Charcoal

```{r}
##parameters [THESE VALUES WERE EXTRACTED FROM KHAVARI ET AL 2023 (tables 12, 13, and 14 in the supplementary material)]
Ucharcoal = 30 #MJ/Kg
CO2_charcoal = 121 #Kg/GJ
CH4_charcoal = 0.5760 #Kg/GJ
N2O_charcoal = 0.0010 #Kg/GJ
BC_charcoal = 0.1075 #Kg/GJ
OC_charcoal = 0.3080 #Kg/GJ

GWP_CH4_charcoal = 25
GWP_N2O_charcoal = 298
GWP_BC_charcoal = 900
GWP_OC_charcoal = -46

##
gamma_charcoal = CO2_charcoal + CH4_charcoal*GWP_CH4_charcoal + N2O_charcoal*GWP_N2O_charcoal + BC_charcoal*GWP_BC_charcoal + OC_charcoal*GWP_OC_charcoal ## emissions unit is kg/GJ

efficiency_traditional_charcoal = .20
efficiency_ICS = .30

##
energyuse_traditional_charcoal = (3.64/efficiency_traditional_charcoal)*10^-3
energyuse_ICS = (3.64/efficiency_ICS)*10^-3

##
emissions_traditional_charcoal = energyuse_traditional_charcoal*gamma_charcoal
emissions_ICS = energyuse_ICS*gamma_charcoal


###additional emission associated with the production of charcoal for each meal

mass_charcoal_used = energyuse_traditional_charcoal*1000/Ucharcoal ##mass of charcoal used for each meal

CO2_charcoal_manufacture = 1.626 #kg
CH4_charcoal_manufacture = 0.0396
BC_charcoal_manufacture = 0.02*10^-3
OC_charcoal_manufacture = 0.74*10^-3

emission_charcoal_manufacturing = mass_charcoal_used*(CO2_charcoal_manufacture + CH4_charcoal_manufacture *GWP_CH4_charcoal + BC_charcoal_manufacture*GWP_BC_charcoal + OC_charcoal_manufacture*10^-3*GWP_OC_charcoal)

##
emissions_traditional_charcoal_total = energyuse_traditional_charcoal*gamma_charcoal + emission_charcoal_manufacturing
emissions_ICS_total = energyuse_ICS*gamma_charcoal + emission_charcoal_manufacturing


##
wt_emissions_traditional_charcoal_total = 1
wt_emissions_ICS_total = 0
##
mean_emissions_charcoal = wt_emissions_traditional_charcoal_total*emissions_traditional_charcoal_total + emissions_ICS_total*wt_emissions_ICS_total
```


## LPG

```{r}

##parameters [THESE VALUES WERE EXTRACTED FROM KHAVARI ET AL 2023 (tables 12, 13, and 14 in the supplementary material)]
ULPG = 45.5 #MJ/Kg
CO2_LPG = 63 #Kg/GJ
CH4_LPG = 0.0028 #Kg/GJ
N2O_LPG = 0.0001 #Kg/GJ
BC_LPG = 0.0044 #Kg/GJ
OC_LPG = 0.0091 #Kg/GJ

GWP_CH4_LPG = 25
GWP_N2O_LPG = 298
GWP_BC_LPG = 900
GWP_OC_LPG = -46



##
gamma_LPG = CO2_LPG + CH4_LPG*GWP_CH4_LPG + N2O_LPG*GWP_N2O_LPG + BC_LPG*GWP_BC_LPG + OC_LPG*GWP_OC_LPG ## emissions unit is kg/GJ

##
efficiency_LPG= .5

##
energyuse_LPG= (3.64/efficiency_LPG)*10^-3


##emission during LPG transportation


##
emissions_LPG = energyuse_LPG*gamma_LPG
```

##electricity

```{r}


##parameters [THESE VALUES WERE EXTRACTED FROM KHAVARI ET AL 2023 (tables 12, 13, and 14 in the supplementary material)]
Uelectricity = 80 #MJ/Kwh

#electricity mix (rpoduced in PJ)
coal_elec = 0.76
oil_elec = 19.87
hydro_elec = 3.82
wind_elec = 0.54
solar_elec =  0.06
  
###
emission_coal_elec= 90 #Kton CO2eq/PJ
emission_oil_elec = 74.69
emission_hydro_elec = 0
emission_wind_elec = 0
emission_solar_elec = 0


##
gamma_electricity = (emission_coal_elec*coal_elec + oil_elec*emission_oil_elec )/(coal_elec + oil_elec + hydro_elec + wind_elec + solar_elec)## emissions unit is kg/PJ

#convert to Kg/MJ

gamma_electricity = gamma_electricity/10^3  #conversion taken for accounting for kton (mass) and PJ(energy)

##
efficiency_electricity= .8

##
energyuse_electricity= (3.64/efficiency_electricity)


##emission during electricity transportation


##
emissions_electricity = energyuse_electricity*gamma_electricity

```


## Data on poulation and urban-rural divide and convert it to a dataframe

```{r}

geoTiffPath_pop <- "/Users/yusufjameel/Dropbox/ProjectDrawdown/BlackCarbon/Open-Source-Spatial-Clean-Cooking-Tool-OnStove-1a72a05/example/SEN1/Demographics/Population/Population.tif"
Populationraster <- raster(geoTiffPath_pop)
Populationraster_df <- as.data.frame(Populationraster,  xy = TRUE)




geoTiffPath_urban_rural <- "/Users/yusufjameel/Dropbox/ProjectDrawdown/BlackCarbon/Open-Source-Spatial-Clean-Cooking-Tool-OnStove-1a72a05/example/SEN1/Demographics/Urban_rural_divide/Urban_rural_divide.tif"
urban_ruralraster <- raster(geoTiffPath_urban_rural)
urban_ruralraster_df <- as.data.frame(urban_ruralraster,  xy = TRUE)

urban_ruralraster_df <- urban_ruralraster_df %>%
  mutate(urban = case_when(
    Urban_rural_divide %in% c(10, 11, 12, 13, 21) ~ 1,
    Urban_rural_divide %in% c(22, 23, 30) ~ 2,
    TRUE ~ NA_integer_  # Optional, assigns NA to other cases
  ))

pop_urban_rural = merge(Populationraster_df, urban_ruralraster_df, by.x = c("x","y"), by.y =c("x", "y"))

```


## add mean emissions 

```{r}

### rural and urban mean emission per meal [all numbers obtained from Khavari et al 2023]

biomass_proportion_rural = 0.81 #### 81% of the fuel used in rural area is biomass
LPG_proportion_rural = 0.03 #
charcoal_proportion_rural = 0.15

biomass_proportion_urban = 0.18
LPG_proportion_urban = 0.34
charcoal_proportion_urban = 0.45
electricity_proportion_urban = 0.01

mean_emission_rural = biomass_proportion_rural*mean_emissions_biomass + LPG_proportion_rural*emissions_LPG + charcoal_proportion_rural*mean_emissions_charcoal

mean_emission_urban = biomass_proportion_urban*mean_emissions_biomass + LPG_proportion_urban*emissions_LPG + charcoal_proportion_urban *mean_emissions_charcoal + electricity_proportion_urban*emissions_electricity


meals_per_day =3
days_year=365


#per meal_emission_per year per fuel type
emissions_biomass_year_household = mean_emissions_biomass*meals_per_day*days_year
emissions_LPG_year_household = emissions_LPG*meals_per_day*days_year
emissions_charcoal_year_household = mean_emissions_charcoal*meals_per_day*days_year
emissions_electricity_year_household = emissions_electricity*meals_per_day*days_year

mean_emission_rural_year_household  = mean_emission_rural*meals_per_day*days_year
mean_emission_urban_year_household  = mean_emission_urban*meals_per_day*days_year

```


## emissions per unit raster. 

##**PLEASE NOTE in the following code **##

##urban == 1 means rural area. 12 is average household rural household size
##urban == 2 means urban area. 11.3 is average household urban household size
## emission values are typically divided by 10^9 to convert the values from kg to million tons

```{r}
pop_urban_rural <- pop_urban_rural %>%
  mutate(mean_emission = case_when(
    urban == 1 ~ mean_emission_rural_year_household*Population/12, 
    urban == 2 ~ mean_emission_urban_year_household*Population/11.3,  
    TRUE ~ NA_real_  # If neither condition is met, assign NA
  ))

total_unclean_emission = sum(pop_urban_rural$mean_emission, na.rm=T)/10^9


####let's look at biomass emissions

pop_urban_rural <- pop_urban_rural %>%
  mutate(mean_emission_biomass = case_when(
    urban == 1 ~ mean_emission_rural_year_household*Population*biomass_proportion_rural/12,
    urban == 2 ~ mean_emission_urban_year_household*Population*biomass_proportion_urban/11.3,
    TRUE ~ NA_real_  # If neither condition is met, assign NA
  ))

###total emissions (urabn and rural)
summary_mean_emission <- pop_urban_rural %>%
  group_by(urban) %>%
  summarize(mean_mean_emission = sum(mean_emission, na.rm = TRUE))

##total emissions [rural and urban]
summary_mean_emission$mean_mean_emission = summary_mean_emission$mean_mean_emission/10^9
# View the summary
print(summary_mean_emission)

summary_pop <- pop_urban_rural %>%
  group_by(urban) %>%
  summarize(summary_pop = sum(Population, na.rm = TRUE))

```

###total emission by fuel type for urban and rural regions [based on population data from the raster]

##**PLEASE NOTE in the following code **##

## 12 is average household rural household size
## 11.3 is average household urban household size
## emission values are typically divided by 10^9 to convert the values from kg to million tons

```{r}
total_emission_biomass_rural = summary_pop$summary_pop[1]*biomass_proportion_rural*emissions_biomass_year_household/10^9/12
total_emission_charcoal_rural = summary_pop$summary_pop[1]*charcoal_proportion_rural*emissions_charcoal_year_household/10^9/12
total_emission_LPG_rural = summary_pop$summary_pop[1]*LPG_proportion_rural*emissions_LPG_year_household/10^9/12

total_emission_biomass_urban = summary_pop$summary_pop[2]*biomass_proportion_urban*emissions_biomass_year_household/10^9/11.3
total_emission_charcoal_urban = summary_pop$summary_pop[2]*charcoal_proportion_urban*emissions_charcoal_year_household/10^9/11.3
total_emission_LPG_urban = summary_pop$summary_pop[2]*LPG_proportion_urban*emissions_LPG_year_household/10^9/11.3
total_emission_electricity_urban = summary_pop$summary_pop[2]*electricity_proportion_urban*emissions_electricity_year_household/10^9/11.3
```



###total emission by fuel type for urban and rural regions [based on population data from the socio-techno file]
##**PLEASE NOTE in the following code **##

## 12 is average household rural household size
## 11.3 is average household urban household size
## emission values are typically divided by 10^9 to convert the values from kg to million tons

```{r}

pop = 16744000
urban = 0.48
rural = 0.52

pop_urban = pop*urban
pop_rural = pop*rural

total_emission_biomass_rural_pop2 = pop_rural*biomass_proportion_rural*emissions_biomass_year_household/10^9/12

total_emission_charcoal_rural_pop2 = pop_rural*charcoal_proportion_rural*emissions_charcoal_year_household/10^9/12

total_emission_LPG_rural_pop2 = pop_rural*LPG_proportion_rural*emissions_LPG_year_household/10^9/12


total_emission_biomass_urban_pop2 = pop_urban*biomass_proportion_urban*emissions_biomass_year_household/10^9/11.3

total_emission_charcoal_urban_pop2 = pop_urban*charcoal_proportion_urban*emissions_charcoal_year_household/10^9/11.3

total_emission_LPG_urban_pop2 = pop_urban*LPG_proportion_urban*emissions_LPG_year_household/10^9/11.3

total_emission_electricity_urban_pop2 = pop_urban*electricity_proportion_urban*emissions_electricity_year_household/10^9/11.3

rural_emission_pop2 = sum(total_emission_biomass_rural_pop2 + total_emission_charcoal_rural_pop2 +total_emission_LPG_rural_pop2 )

urban_emision_pop2 = sum(total_emission_biomass_urban_pop2+total_emission_charcoal_urban_pop2+total_emission_LPG_urban_pop2+total_emission_electricity_urban_pop2)
```

## emissions per unit raster if all switched to electric or all LPG

```{r}

#electric
pop_urban_rural <- pop_urban_rural %>%
  mutate(all_electric = case_when(
    urban == 1 ~ Population*meals_per_day*days_year*emissions_electricity/12,
    urban == 2 ~ Population*meals_per_day*days_year*emissions_electricity/11.3,
    TRUE ~ NA_real_  # If neither condition is met, assign NA
  ))


emissions_all_electric = sum(pop_urban_rural$all_electric ,na.rm=T)/10^9

##LPG
pop_urban_rural <- pop_urban_rural %>%
  mutate(all_LPG = case_when(
    urban == 1 ~ Population*meals_per_day*days_year*emissions_LPG/12,
    urban == 2 ~ Population*meals_per_day*days_year*emissions_LPG/11.3,
    TRUE ~ NA_real_  # If neither condition is met, assign NA
  ))


emissions_all_LPG = sum(pop_urban_rural$all_LPG ,na.rm=T)/10^9
```

##final emissions based on Khavari estimates (quick estimate)

```{r}

## total population of senegal is 16564068 [in the socioeconomic file it is 16744000]
# It assumes 72% of the population uses electricity
# It assumes 28% of the population uses LPG

emissions_senegal_clean = ((pop_urban*0.72*emissions_electricity_year_household/11.3 + pop_rural*0.72*emissions_electricity_year_household/12) + (pop_urban*0.28*emissions_LPG_year_household/11.3 + pop_rural*0.28*emissions_LPG_year_household/12))/10^9
emissions_senegal_clean
```

## plots

#plot 1. emissions per meal by fuel category

```{r}


data <- data.frame(
  Variable = c("Traditional Biomass", "Charcoal", "LPG", "Electricity"),
  Value = c(mean_emissions_biomass, mean_emissions_charcoal, emissions_LPG, emissions_electricity)
)

data$Variable <- factor(data$Variable, levels = data$Variable[order(data$Value)])

# Create a histogram
ggplot(data, aes(x = Variable, y = Value)) +
  geom_bar(stat = "identity", fill = "skyblue", color = "black") +
  labs(x = "Fuel Type", y = "Kg CO2-eq emitted per meal", title = "GHG emissions from cooking fuel for Senegal") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.border = element_rect(color = "black", fill = NA, size = 1))  # Rotate x-axis labels for better readability

```

##plot 2 total emissions from each fuel type 

```{r}

data1 <- data.frame(
  Region = rep("Rural", 3),
  Category = c("Biomass", "Charcoal", "LPG"),
  Proportion = c(biomass_proportion_rural,charcoal_proportion_rural, LPG_proportion_rural),
  Emission = c(total_emission_biomass_rural, total_emission_charcoal_rural , total_emission_LPG_rural)
)

# Create a dataframe for the second histogram
data2 <- data.frame(
  Region = rep("Urban", 4),
  Category = c("Biomass", "Charcoal", "LPG", "Electricity"),
  Proportion = c(biomass_proportion_urban,charcoal_proportion_urban, LPG_proportion_urban, electricity_proportion_urban),
  Emission = c(total_emission_biomass_urban, total_emission_charcoal_urban, total_emission_LPG_urban, total_emission_electricity_urban))

# Combine the dataframes
combined_data <- rbind(data1, data2)

# Create a histogram
ggplot(combined_data, aes(fill = Category, y = Emission, x = Region)) +
  geom_bar(position = "stack", stat = "identity", color = "black", size = 0.2) +
  labs(x = "Region", y = "Total Emission - Mton CO2eq", title = "Total Emission by fuel type and Region") +
  scale_fill_manual(values = c("Biomass" = "red", "Charcoal" = "yellow", "LPG" = "lightgreen", "Electricity" = "lightblue")) +
  theme_minimal()+
    theme(panel.border = element_rect(color = "black", fill = NA, size = 0.5),
        axis.line = element_line(color = "black"),
        legend.position = "right",
        legend.title = element_blank())

```

##plot3: map of emissions - current emission from Senegal

#step 1: get level1 and level 3 boundaries map

```{r}

## import level 3 and level 1 maps from Senegal

#import level 3 administrative boundaries 

Senegallevel3<- "/Users/yusufjameel/Dropbox/ProjectDrawdown/BlackCarbon/Open-Source-Spatial-Clean-Cooking-Tool-OnStove-1a72a05/example/SEN1/GADM41Rasters/SEN_GADM3_KhaveriGrid.tif"
Senegallevel3<- raster(Senegallevel3)
Senegallevel3_df = as.data.frame(Senegallevel3, xy=T)
Senegallevel3[Senegallevel3 == 0] <- NA
shp_Senegallevel3 <- rasterToPolygons(Senegallevel3, dissolve = T)
shp_Senegallevel3_sf <- sf::st_as_sf(shp_Senegallevel3)

#import level 1 administrative boundaries 

Senegallevel1<- "/Users/yusufjameel/Dropbox/ProjectDrawdown/BlackCarbon/Open-Source-Spatial-Clean-Cooking-Tool-OnStove-1a72a05/example/SEN1/GADM41Rasters/SEN_GADM1_KhaveriGrid.tif"
Senegallevel1<- raster(Senegallevel1)
Senegallevel1_df = as.data.frame(Senegallevel1, xy=T)
Senegallevel1[Senegallevel1 == 0] <- NA
shp_Senegallevel1 <- rasterToPolygons(Senegallevel1, dissolve = T)
shp_Senegallevel1_sf <- sf::st_as_sf(shp_Senegallevel1)


```


##step 2: merge emissions df with GADM level 3 rasters

```{r}
##merge Senegal level 3 and level 1 data sets

merge_level3_1_df <- merge(Senegallevel1_df,Senegallevel3_df, by.x =c("x","y"), by.y =c("x","y"))

###merge level 3 and level 1 data sets with the dataset with all the calculations

pop_urban_rural_spatial = merge(merge_level3_1_df,fuel_senegal, by.x =c("x","y"), by.y =c("x","y"))
```

##step 3: summarize total emissions at level 3 boundaries

```{r}

##summarize total emission at level 3 boundary

emissions_level3 <- pop_urban_rural_spatial  %>%
  group_by(SEN_GADM3_KhaveriGrid) %>%
  summarise(Sum_Value = sum(mean_emission, na.rm=T))


emissions_level3 <- emissions_level3 %>%
  rename(CarbonEmissions = Sum_Value)

emissions_level3$CarbonEmissions = emissions_level3$CarbonEmissions/10^9 #convert kg to million ton
emissions_level3 <- as.data.frame(emissions_level3)


##merge summarized emissions at level 3 with level 3 raster

##merge results to the Senegal level 2 raster
emissionslevel3_raster_df = merge(Senegallevel3_df, emissions_level3, by = "SEN_GADM3_KhaveriGrid", all=T)
emissionslevel3_raster_df  = emissionslevel3_raster_df[emissionslevel3_raster_df$SEN_GADM3_KhaveriGrid != 0,]


##export the emissions at level 3 boundaries
 
emissionslevel3_raster_df_export <- emissionslevel3_raster_df %>%
  select(x, y, CarbonEmissions)
raster_emission_level3 <- raster::rasterFromXYZ(emissionslevel3_raster_df_export)
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #### projection
proj4string(raster_emission_level3) = wgs.84
raster_emission_level3

writeRaster(raster_emission_level3, filename = "reduction_in_GHG_emissions.tif", format = "GTiff")

```


##step 4:  create a color palette and plot total emissions

```{r}
emissionslevel3_raster_df$Category <- cut(emissionslevel3_raster_df$CarbonEmissions, breaks = c(0, 0.004, 0.01, 0.02, 0.04,0.06, 0.08, 0.1, 0.2, 0.3, 0.7), labels = c("<0.004", "0.004-0.01", "0.01-0.02","0.02-0.04", "0.04-0.06","0.06-0.08", "0.08-0.1", "0.1-0.2", "0.2-0.3", "0.3-0.4"))

#color_palette =brewer.pal(n = 9, name = 'YlOrRd')

#category_colors <- c("<0.004" = color_palette[1], "0.004-0.01" = color_palette[2], "0.01-0.02" = color_palette[3], "0.02-0.04" = color_palette[4],"0.04-0.06" = color_palette[5], "0.06-0.08" = color_palette[6],"0.08-0.1" = color_palette[7], "0.1-0.2" = color_palette[8],"0.2-0.4" = color_palette[9] )


category_colors <- c("<0.004" = "gray98","0.004-0.01" = "#e5f3fd", "0.01-0.02"= "#bdd5e7","0.02-0.04"= "#9abddc", "0.04-0.06" = "wheat","0.06-0.08"= "gold","0.08-0.1" = "orange2","0.1-0.2"= "red1","0.2-0.3" ="red3","0.3-0.4" = "maroon4")

 ggplot() +
        geom_raster(data = emissionslevel3_raster_df, aes(x=x, y=y, fill = Category)) +
   scale_fill_manual(values = category_colors, name = "carbon emissions (MT)") +
      #geom_sf(data = shp_senegallevel2_sf, fill = NA, color = "gray75", linetype = "dashed") +
  geom_sf(data = shp_Senegallevel3_sf, fill = NA, color = "gray40", linewidth = 0.2, linetype = "longdash") +
  geom_sf(data = shp_Senegallevel1_sf, fill = NA, color = "black", linewidth = 0.4) +
  theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 0.4)) +
  ggtitle("Annual GHG emission") +
  xlab("Longitude") + ylab("Latitdue")
 
```

##Plot 4: emission from all electric

#step 1: summarize electric emissions at GADM level 3 

```{r}
electric_level3 <- pop_urban_rural_spatial  %>%
  group_by(SEN_GADM3_KhaveriGrid) %>%
  summarise(Sum_Value = sum(all_electric, na.rm=T))


electric_level3 <- electric_level3 %>%
  rename(CarbonEmissions = Sum_Value)

electric_level3$CarbonEmissions = electric_level3$CarbonEmissions/10^9 #convert kg to million ton
electric_level3 <- as.data.frame(electric_level3)


##merge summarized electric at level 3 with level 3 raster

##merge results to the Senegal level 2 raster
electriclevel3_raster_df = merge(Senegallevel3_df, electric_level3, by = "SEN_GADM3_KhaveriGrid", all=T)
electriclevel3_raster_df  = electriclevel3_raster_df[electriclevel3_raster_df$SEN_GADM3_KhaveriGrid != 0,]

##export the emissions at level 3 boundaries
 
electriclevel3_raster_df_export <- electriclevel3_raster_df %>%
  select(x, y, CarbonEmissions)
raster_electric_level3 <- raster::rasterFromXYZ(electriclevel3_raster_df_export)
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #### projection
proj4string(raster_electric_level3) = wgs.84
raster_electric_level3

writeRaster(raster_electric_level3, filename = "GHG_emissions_electric.tif", format = "GTiff")

```

##step 2: create a color palette and plot electric emission

```{r}

electriclevel3_raster_df$Category <- cut(electriclevel3_raster_df$CarbonEmissions, breaks = c(0, 0.004, 0.01, 0.02, 0.04,0.06, 0.08, 0.1, 0.2, 0.3, 0.7), labels = c("<0.004", "0.004-0.01", "0.01-0.02","0.02-0.04", "0.04-0.06","0.06-0.08", "0.08-0.1", "0.1-0.2", "0.2-0.3", "0.3-0.4"))

#color_palette =brewer.pal(n = 9, name = 'YlOrRd')

#category_colors <- c("<0.004" = color_palette[1], "0.004-0.01" = color_palette[2], "0.01-0.02" = color_palette[3], "0.02-0.04" = color_palette[4],"0.04-0.06" = color_palette[5], "0.06-0.08" = color_palette[6],"0.08-0.1" = color_palette[7], "0.1-0.2" = color_palette[8],"0.2-0.4" = color_palette[9] )


category_colors <- c("<0.004" = "gray98","0.004-0.01" = "#e5f3fd", "0.01-0.02"= "#bdd5e7","0.02-0.04"= "#9abddc", "0.04-0.06" = "wheat","0.06-0.08"= "gold","0.08-0.1" = "orange2","0.1-0.2"= "red1","0.2-0.3" ="red3","0.3-0.4" = "maroon4")

 ggplot() +
        geom_raster(data = electriclevel3_raster_df, aes(x=x, y=y, fill = Category)) +
   scale_fill_manual(values = category_colors, name = "carbon emissions (MT)") +
      #geom_sf(data = shp_senegallevel2_sf, fill = NA, color = "gray75", linetype = "dashed") +
  geom_sf(data = shp_Senegallevel3_sf, fill = NA, color = "gray40", linewidth = 0.2, linetype = "longdash") +
  geom_sf(data = shp_Senegallevel1_sf, fill = NA, color = "black", linewidth = 0.4) +
  theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 0.4)) +
  ggtitle("Annual GHG emission - all electric source") +
  xlab("Longitude") + ylab("Latitdue")

```

##Plot 5: emission from all LPG

#step 1: summarize electric emissions at GADM level 3 

```{r}
LPG_level3 <- pop_urban_rural_spatial  %>%
  group_by(SEN_GADM3_KhaveriGrid) %>%
  summarise(Sum_Value = sum(all_LPG, na.rm=T))


LPG_level3 <- LPG_level3 %>%
  rename(CarbonEmissions = Sum_Value)

LPG_level3$CarbonEmissions = LPG_level3$CarbonEmissions/10^9 #convert kg to million ton
LPG_level3 <- as.data.frame(LPG_level3)


##merge summarized LPG at level 3 with level 3 raster

##merge results to the Senegal level 2 raster
LPGlevel3_raster_df = merge(Senegallevel3_df, LPG_level3, by = "SEN_GADM3_KhaveriGrid", all=T)
LPGlevel3_raster_df  = LPGlevel3_raster_df[LPGlevel3_raster_df$SEN_GADM3_KhaveriGrid != 0,]

##export the emissions at level 3 boundaries
 
LPGlevel3_raster_df_export <- LPGlevel3_raster_df %>%
  select(x, y, CarbonEmissions)
raster_LPG_level3 <- raster::rasterFromXYZ(LPGlevel3_raster_df_export)
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #### projection
proj4string(raster_LPG_level3) = wgs.84
raster_LPG_level3

writeRaster(raster_LPG_level3, filename = "GHG_emissions_LPG.tif", format = "GTiff")


```

##step 2: create a color palette and plot LPG emission

```{r}

LPGlevel3_raster_df$Category <- cut(LPGlevel3_raster_df$CarbonEmissions, breaks = c(0, 0.004, 0.01, 0.02, 0.04,0.06, 0.08, 0.1, 0.2, 0.3, 0.7), labels = c("<0.004", "0.004-0.01", "0.01-0.02","0.02-0.04", "0.04-0.06","0.06-0.08", "0.08-0.1", "0.1-0.2", "0.2-0.3", "0.3-0.4"))

#color_palette =brewer.pal(n = 9, name = 'YlOrRd')

#category_colors <- c("<0.004" = color_palette[1], "0.004-0.01" = color_palette[2], "0.01-0.02" = color_palette[3], "0.02-0.04" = color_palette[4],"0.04-0.06" = color_palette[5], "0.06-0.08" = color_palette[6],"0.08-0.1" = color_palette[7], "0.1-0.2" = color_palette[8],"0.2-0.4" = color_palette[9] )


category_colors <- c("<0.004" = "gray98","0.004-0.01" = "#e5f3fd", "0.01-0.02"= "#bdd5e7","0.02-0.04"= "#9abddc", "0.04-0.06" = "wheat","0.06-0.08"= "gold","0.08-0.1" = "orange2","0.1-0.2"= "red1","0.2-0.3" ="red3","0.3-0.4" = "maroon4")

 ggplot() +
        geom_raster(data = LPGlevel3_raster_df, aes(x=x, y=y, fill = Category)) +
   scale_fill_manual(values = category_colors, name = "carbon emissions (MT)") +
      #geom_sf(data = shp_senegallevel2_sf, fill = NA, color = "gray75", linetype = "dashed") +
  geom_sf(data = shp_Senegallevel3_sf, fill = NA, color = "gray40", linewidth = 0.2, linetype = "longdash") +
  geom_sf(data = shp_Senegallevel1_sf, fill = NA, color = "black", linewidth = 0.4) +
  theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 0.4)) +
  ggtitle("Annual GHG emission - all LPG source") +
  xlab("Longitude") + ylab("Latitdue")

```

##plot 5: Histogram of emissions under different scenarios

```{r}

total_current_emission = sum(combined_data$Emission)
total_emission_electricity = emissions_all_electric
total_emission_LPG = emissions_all_LPG
total_emission_max_benefit = emissions_max_benefit



data <- data.frame(
  Variable = c("Current fuel mix", "Electricity", "LPG", "Mix with maximum benefit"),
  Value = c(total_current_emission, total_emission_electricity, total_emission_LPG, total_emission_max_benefit)
)

# Create a histogram
ggplot(data, aes(x = Variable, y = Value)) +
  geom_bar(stat = "identity", fill = "skyblue", color = "black") +
  labs(x = "Fuel Type", y = "Kg CO2-eq emitted per year", title = "GHG emissions from different fuel mix in Senegal") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.border = element_rect(color = "black", fill = NA, size = 1))  # Rotate x-axis labels for better readability


``` 

##plot 6: Lives saved from implementing the fuel mix with maximum benefit

```{r}

MortalityTiffPath <- "/Users/yusufjameel/Dropbox/ProjectDrawdown/BlackCarbon/Open-Source-Spatial-Clean-Cooking-Tool-OnStove-1a72a05/example/SEN1/Rasters/deaths_avoided_per_100k.tif"
Mortalityraster <- raster(MortalityTiffPath)

Mortality_df <- as.data.frame(Mortalityraster,  xy = TRUE)

Mortality_df_Senegal_pop = merge(Mortality_df, pop_urban_rural_spatial, by.x =c("x","y"), by.y =c("x","y"))


###get mean mortality rate  per level 3 district
result_mortality <- Mortality_df_Senegal_pop  %>%
  group_by(SEN_GADM3_KhaveriGrid) %>%
  summarise(mean = mean(deaths_avoided_per_100k, na.rm=T))

result_mortality <- result_mortality %>%
  rename(Deathsavoided = mean)


result_pop <- Mortality_df_Senegal_pop   %>%
 group_by(SEN_GADM3_KhaveriGrid) %>%
  summarise(sum_Value = sum(Population, na.rm=T))


##get the population per level 3 district
result_pop <- result_pop %>%
  rename(Population = sum_Value)

result_pop <- as.data.frame(result_pop)

result_total_mortality = merge(result_pop, result_mortality, by = "SEN_GADM3_KhaveriGrid", all=T)

###this calculates the total number of avoided deaths per district
result_total_mortality$totalmortality =  result_total_mortality$Population*result_total_mortality$Deathsavoided/100000


####merge with the raster


totalmortalitylevel3 = merge(Senegallevel3_df, result_total_mortality, by = "SEN_GADM3_KhaveriGrid", all=T)
totalmortalitylevel3  = totalmortalitylevel3 [totalmortalitylevel3$SEN_GADM3_KhaveriGrid != 0,]

##export the emissions at level 3 boundaries
 
totalmortalitylevel3_df_export <- totalmortalitylevel3 %>%
  select(x, y, totalmortality)
totalmortalitylevel3_df_export  <- raster::rasterFromXYZ(totalmortalitylevel3_df_export)
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #### projection
proj4string(totalmortalitylevel3_df_export) = wgs.84
totalmortalitylevel3_df_export

#writeRaster(totalmortalitylevel3_df_export, filename = "Avoided_death.tif", format = "GTiff")


###

totalmortalitylevel3$Category <- cut(totalmortalitylevel3$totalmortality, breaks = c(0, 20, 40, 60, 80,100, 150, 200, 400, 500), labels = c("<20", "20-40", "40-60","60-80", "80-100","100-150", "150-200", "200-400", "400-500"))


color_palette =brewer.pal(n = 9, name = 'Blues')

category_colors <- c("<20" = color_palette[1], "20-40" = color_palette[2], "40-60" = color_palette[3], "60-80" = color_palette[4],"80-100" = color_palette[5], "100-150" = color_palette[6],"150-200" = color_palette[7], "200-400" = color_palette[8],"400-500" = color_palette[9] )


# Create a ggplot object and plot the specified field
ggplot() +
        geom_raster(data = totalmortalitylevel3, aes(x=x, y=y, fill = Category)) +
   scale_fill_manual(values = category_colors, name = "Lives saved") +
  geom_sf(data = shp_Senegallevel3_sf, fill = NA, color = "gray40", linewidth = 0.2, linetype = "longdash") +
  geom_sf(data = shp_Senegallevel1_sf, fill = NA, color = "black", linewidth = 0.4) +
  theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 0.4)) +
  ggtitle("Lives saved after switching to cleaner fuels") +
  xlab("Longitude") + ylab("Latitdue")

```


## plot 7: time saved

```{r}
TimesavedTiffPath <- "/Users/yusufjameel/Dropbox/ProjectDrawdown/BlackCarbon/Open-Source-Spatial-Clean-Cooking-Tool-OnStove-1a72a05/example/SEN1/Rasters/time_saved_per_household.tif"
Timesavedraster <- raster(TimesavedTiffPath)


Timesavedraster_df <- as.data.frame(Timesavedraster ,  xy = TRUE)


Timesavedraster_df_Senegal = merge(Timesavedraster_df, pop_urban_rural_spatial, by.x =c("x","y"), by.y =c("x","y"))

result_timesaved <- Timesavedraster_df_Senegal  %>%
  group_by(SEN_GADM3_KhaveriGrid) %>%
  summarise(mean = mean(time_saved_per_household, na.rm=T))

result_timesaved <- result_timesaved %>%
  rename(Timesaved = mean)


result_timesaved <- as.data.frame(result_timesaved)



##merge results to the Senegal level 3 raster
timesavedlevel3 = merge(Senegallevel3_df, result_timesaved, by = "SEN_GADM3_KhaveriGrid", all=T)
timesavedlevel3   = timesavedlevel3[timesavedlevel3$SEN_GADM3_KhaveriGrid != 0,]

##########

timesavedlevel3_df_export <- timesavedlevel3 %>%
  select(x, y, Timesaved)
timesavedlevel3_df_export  <- raster::rasterFromXYZ(timesavedlevel3_df_export)
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #### projection
proj4string(timesavedlevel3_df_export) = wgs.84
timesavedlevel3_df_export

#writeRaster(timesavedlevel3_df_export, filename = "Time_saved.tif", format = "GTiff")


########
timesavedlevel3$Category <- cut(timesavedlevel3$Timesaved, breaks = c(0, .20, .40, .60, .80,1.00, 1.20, 1.400, 1.6, 1.8), labels = c("<0.20", "0.20-0.40", "0.4-0.6","0.6-0.8", "0.8-1.0","1.0-1.2", "1.2-1.4", "1.4-1.6", "1.6-1.8"))


color_palette =brewer.pal(n = 9, name = 'Blues')

category_colors <- c("<0.20" = color_palette[1], "0.20-0.40" = color_palette[2], "0.4-0.6" = color_palette[3], "0.6-0.8" = color_palette[4],"0.8-1.0" = color_palette[5], "1.0-1.2" = color_palette[6],"1.2-1.4" = color_palette[7], "1.4-1.6" = color_palette[8],"1.6-1.8" = color_palette[9] )





# Create a ggplot object and plot the specified field
ggplot() +
        geom_raster(data = timesavedlevel3, aes(x=x, y=y, fill = Category)) +
   scale_fill_manual(values = category_colors, name = "Time saved (hours)") +
  geom_sf(data = shp_Senegallevel3_sf, fill = NA, color = "gray40", linewidth = 0.2, linetype = "longdash") +
  geom_sf(data = shp_Senegallevel1_sf, fill = NA, color = "black", linewidth = 0.4) +
  theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 0.4)) +
  ggtitle("Time saved per household after switching to cleaner fuels") +
  xlab("Longitude") + ylab("Latitdue")

```


